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FOREWORD 


Opinions, interpretations, conclusions and recommendations are 
those of the author and are not necessarily endorsed by the U.S. 
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material. 
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this report do not constitute an official Department of the Army 
endorsement or approval of the products or services of these 
organizations. 
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adhered to the "Guide for the Care and Use of Laboratory Animals," 
prepared by the Committee on Care and Use of Laboratory Animals of 
the Institute of Laboratory Animal Resources, National Research 
Council (NIH Publication No. 86-23, Revised 1985). 


( ) For the protection of human subjects, the investigator(s) 
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Introduction 

The military employs a broad range of laser radiation in training devices, rangefinders, 
target designators, communications devices and other instruments. This equipment emits either 
single pulses or sequences of pulses in beams of various diameters. The research performed under 
this contract directly supports the U. S. Army Medical Research and Materiel Command 
(USAMRMC) mission to assess the health effects and hazards of nonionizing electromagnetic 
radiation from such laser systems. The data obtained will support evaluation of current permissible 
exposure limits and will aid health policy makers, both within and outside the DoD, in developing 
injury prevention criteria. The general approach is to make direct determinations of damage 
threshold levels for non-ionizing radiation for specific exposure conditions (e.g., wavelengths, 
pulse durations, etc.) and to develop models of the damage mechanism that enable the extension of 
the results to other exposure conditions. 

Under past support from the Army Medical Research and Development Command we 
determined corneal damage thresholds for CO 2 lasers that emit single and sequences of pulses 
having durations of 1 ms and longer and developed thermal damage models for predicting this type 
of damagei-5. We also determined comeal damage thresholds for single 80 ns pulses of CO 2 laser 

radiation and for three multiple-pulse (2 pulses at 1 Hz.; and 2 and 8 pulses at 10 Hz.).®’ ^ In the 
case of these very short pulses, damage mechanisms other than thermal (e.g., acousticd pressure 
pulses) could also play a role. Light and electron microscopy revealed unusual dismptions of the 
anterior epithelial surface for the threshold single-pulse exposure. The characteristics of these 
dismptions differed from those observed with simple thermal damage at longer pulse durations and 
appeared to be consistent with a mechanical [e.g., acoustic] damage mechanism. However, the 
calculated temperature increase produced by the threshold exposure was only slightly lower than 

that calculated for threshold exposures having durations of 1 ms and longer.®- ^ Thus we could not 
exclude a thermal damage mechanism, with the sharp temperature gradients leading to marked 
differences in the character of the damage as compared to damage from longer duration exposures. 

This study aims to expand the database of thresholds for sequences of 80 ns of CO 2 laser 
pulses to include larger numbers of pulses and other pulse repetition frequencies. In addition it 
will attempt to clarify the damage mechanism(s) for these types of pulses and ultimately will begin 
to address damage from mid-infrared wavelengths where the radiation is more penetrating. 

Methods 

Short-pulse exposures are made with a Boston Laser (Model 220S) CO 2 -TEA laser 
operated in the TEMqo mode. This laser delivers 80 ns pulses at pulse repetition frequencies up to 
16 Hz. Mode quality is verified and the beam diameter is measured at the beginning, middle and 
end of each experimental session using a Spiricon linear pyroelectric array. The detector has 64 
elements on 200 |i.m centers. It is mounted on a vertical micropositioner and is read out with a 
LeCroy 9354M digital oscilloscope. Pulse energy is measured with a Scientech detector 
immediately before and immediately after each exposure. 

New Zealand white rabbits of either sex weighing 4-5 pounds are used for the 
experiments. The rabbits are treated in accordance with the Guide for the Care and Use of 
Laboratory Animals (DHEW Publication No. (NIH) 85-23, Revised Edition, 1985 and with the 
Association for Research in Vision and Ophthalmology Resolution on the Use of Animals in 
Research. Prior to exposure the rabbits are anesthetized with an intramuscular injection of xylazine 
and ketamine hydrochloride (Rompun-Ketaset) in the proportions: 60% of 20 m^ml Rompun to 
40% of 100 mg/ml Ketaset by volume. In addition, a topical anesthesia (proparacaine 
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hydrochloride 1/2%) is applied to each eye before exposure and a drop of homatropine bromide 
5% is instilled to dilate the pupil. This facilitates examination of the exposed corneas for minimal 
lesions. The anesthetized animals are placed in a conventional holder for exposure. A speculum is 
inserted in the eye about one minute before exposure and the eye irrigated with BSS solution 
(Alcon Surgical) which is at room temperature; however, in order to create a reproducible “tear 
film,” the irrigation is stopped 20 sec before the exposure and excess fluid is blotted at the limbus. 
The cornea surface is assumed to have returned to its normal temperature after this time. One 
exposure is made to each eye. One-half hour after exposure the rabbits, still under anesthesia, are 
sacrificed with Beuthanasia-D administered in an ear vein. The eyes are enucleated and examined 
for damage using a Nikon FS-3 photo slit-lamp. In selected cases the globes are placed in 
gluteraldehyde/formaldehyde fixative and are delivered to Prof. W. R. Green’s laboratory at the 
Wilmer Ophthalmological Institute where they are processed for microscopy and the micrographs 
evaluated by Dr. Green. 

In a few experiments we exposed enucleated eyes which were at room temperature to test 
the damage mechanism. We have shown previously that reliable damage thresholds can be 

determined in freshly enucleated eyes.^ For these experiments the rabbits were anesthetized and 
given a topical anesthesia and the pupils dilated exactly as was done for the in vivo exposures. The 
rabbits were then sacrificed and their eyes enucleated. The enucleated eyes were placed in room 
temperature (usually ~20 °C) BSS and allowed to equilibrate for at least 5 min. They were then 
exposed using the same protocol as for the in vivo exposures. After exposure the eyes were placed 
back in the BSS solution for 1/2 hr before examining them for damage. 

The criterion that we use for minimal epithelial damage is that due to Brownell and Stuck^, 
namely the presence of a superficial gray-white spot that develops within 1/2 hr after exposure. 
We have found that the damage threshold is sharply defined; i.e., only rarely is there overlap 
between exposures that produce minimal lesions and those that do not. Therefore we do not use 
statistical procedures such as probit analysis in order to determine the threshold, as these would 
require the use of more animals than we deem necessary. We make one exposure per eye, 
bracketing exposures above and below threshold. The bracket is narrowed until there is only about 
a 10% difference in energy between an exposure that produces a minimal lesion and one that does 
not. The threshold exposure is taken to be at the center of the bracket. 

Temperature calculations are based on a Green’s function solution to the heat conduction 
equation for an incident beam with a Gaussian irradiance profile that is absorbed according to 
Beer’s law. The beam is assumed to impinge on a semi-infinite slab and to have a constant peak 
irradiance for the duration of the exposure. We also assume that conduction is the only mode of 
heat transfer and that no heat is lost to the air at the epithelial interface. The absorption coefficient 

and thermal properties are assumed to be those of water.f 2 The solution T{r,z,t), where r is the 
radial distance from the beam axis, z is the depth into the cornea, and t is time, has the form of a 
definite integral that can be evaluated numerically.!- 4, lO since our last Annual report we have 
made additional modifications to this program to facilitate the calculations for multiple-pulse 
exposures using very short pulses. The modifications are discussed below and a copy of the new 
program is listed in Appendix 2. 
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Results and Discussion 

We noted in the first Annual Report covering the period 1 Nov. 1995-31 October 1996 that 
the experimental apparatus, which had not been in use since our previous contact expired in 1989 
and which had been disassembled, was operational and that its operation had been verified on a 
few test exposures. As a result, most of the objectives in the original Statement of Work for the 
first year were shifted to this year. These objectives were: 

1. The data base of threshold conditions for 80 ns pulses from CO 2 -TEA lasers will be 
extended by measuring the threshold energy densities for sequences of 32, 128 and 
1024 pulses at 10 Hz and sequences of 2, 8, 32 and 128 pulses at 20 Hz. The 
existing thresholds for 2 and 8 pulses at 10 Hz will be refined. 

2. An understanding of the damage mechanism for such pulses will be pursued by: 

a) determining how lowering the temperature of the epithelium affects the threshold, 

b) obtaining light and electron micrographs of the damage, 

c) using high-speed photography to investigate ablation of material from the comeal 
surface, and 

d) measuring pressure transients in model systems. 

3. The required characteristics of potential near-infrared laser source(s) for damage 
studies will be identified by performing detailed temperature calculations in the 1.3 to 
2.5 fxm wavelength region. 

This year we have made significant progress toward achieving these objectives in spite of 
having had to move and re-establish our laboratory during the third quarter. The move to 
improved space was necessary because of a reorganization at the Applied Physics Laboratory. 

We repeated and refined the determination of the damage threshold for 2 pulses at 10 Hz. 
The new threshold of 235 mJ/cmVpulse is slightly higher than the value of 200 mJ/cmVpulse that 

we had determined in preliminary experiments.^’ We then confirmed our earlier value for the 
threshold for 8 pulses at 10 Hz. After these experiments we determined epithelial damage 
thresholds for sequences of 32,128, and 1024 pulses at a pulse repetition frequency of 10 Hz. In 
order to perform the experiments with 32 or more pulses it was necessary to put a partMy 
reflecting attenuator in the laser beam because the required pulse energies were below the lasing 
threshold. We found that the beam diameter measured with the partial reflector in place tended to 
be smaller than when it was removed. This led us to question the previously reported single-pulse 

threshold,^’ because at the time it was determined our practice was to attenuate the beam so that 
the Spiracon detector used to measure beam diameter would not be saturated. We noted that if the 
unattenuated beam used in determining the threshold actually had a larger diameter, then the correct 
value of EDu, would be lower. Consequently we repeated the single-pulse threshold measurement 
without using the reflector. In these new experiments the beam intensity was reduced by 
decreasing the laser power. The newly determined threshold is 307 mJ/cm^ compared to the value 
360 mJ/cm^ found previously. The final 10 Hz damage thresholds are listed in Table I together 
with the calculated peak temperature increases. The temperature calculations and their implications 
are discussed below. These threshold results were presented at the 1997 International Laser Safety 
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Conference! 1 and at the annual meeting of the Association for Research in Vision and 
Ophthalmology.! ^ 

In the original Statement of Work we also were to determine epithelial damage thresholds at 
a pulse repetition frequency of 20 Hz. However, we discovered that output of the Boston Laser 
(Model 220S) was unstable at pulse frequencies above 16 Hz (although the specifications claimed 
20 Hz). Boston Laser is no longer in business, consequently technical support for the laser is not 
available. We discussed this problem with our Contracting Officer’s Representative, Mr. Bruce 
Stuck, and it was agreed that we would make the determinations at 16 Hz. The final results are 
listed in Table I (note that we made an additional threshold determination at 1024 pulses which was 
not included in the original Statement of Work). 

Table 1: Threshold Energy Densities and Calcuiated Maximum 
Temperature Rises for Sequences of 80 ns Pulses. 







1 

- 

307 

3.72 

30.25 

2 

10 

235 

3.48 

25.68 

8 

10 

228 

3.80 

32.00 

32 

10 

154 

3.78 

29.15 

128 

10 

136 

3.41 

32.45 

1024 

10 

95 

3.21 

26.60 

2 

16 

265 

3.62 

29.73 

8 

16 

205 

3.75 

31.26 

32 

16 

150 

3.73 

32.90 

128 

16 

105 

3.82 

32.80 

1024 

16 

85 

3.58 

35.06 


*Calculated on the beam axis 10 gm beneath the anterior tear surface. 

The epithelial damage thresholds listed in Table 1 are plotted in Figure 1. The least squares 
fits to these data show that the thresholds at 10 Hz and 16 Hz are correlated by an empirical power 
law of the form 


ED,, = CN-\ 

in which N is the number of pulses in the sequence. The coefficient C and exponent a appear to 
differ slightly for the two cases; however it is not possible to discern if the difference is real. For 
the 10 Hz threshold, C = 291 mJ/cmVpulse and a = 0.162 (R = 0.976), and for the 16 Hz 
threshold, C = 300 mJ/cmVpulse and a = 0.194 (R = 0.997). The least squares fits fall within the 
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±10 percent accuracy estimated from the bracketing procedure used in determining the thresholds. 
If both sets of data are assumed to be part of the same population we find C = 295.5 mJ/cmVpulse 
and a = 0.178 (R = 0.984). The power law is of the same form that we found previously for 
sequences of pulses having individual pulse durations between 0.001 and 1 sec and pulse 
repetition frequencies between 1 and 100 Hz, except that, for these longer pulses, the exponent a 
= 0.25 and the coefficient depended on the duration of the individual pulses.^ Coincidentally, 
retinal damage thresholds for sequences of pulses are also described by a power law of this same 
form with a = 0.25.13 

In order to calculate the temperature histories for these pulse sequences it was necessary to 
make minor modifications to our thermal program in order to obtain sufficient time resolution ^r 
the 80 ns pulses to capture the maximum temperature increase and yet not have too many points for 
a practical calculation. Thus the modification allowed for calculating temperatures at a number of 



Figure 1. The dependence of the threshold energy density per pulse on 
the number of pulses at pulse frequencies of 10 and 16 Hz. The lines are 
least squares fits to a power law of the form ED,^ = CN'". The 

corresponding values of C and a are given in the text. The error bars are 
±10 percent of the values and represent the estimated accuracy of the 
bracketing procedure used in determining the thresholds. 

equally spaced points during each pulse and at a number of points with a wider interval between 
them after each pulse. The numbers of points during and after each pulse and the time interval after 
each pulse for which temperatures are calculated are inputs into the revised program. A listing of 
the revised program is given in Appendix 2. 

We have completed the temperature calculations at the damage threshold for all of the 10 Hz 
and 16 Hz exposures. Figure 2 shows typical temperature histories at the damage threshold for 
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sequences of 32 and 128 pulses at a pulse repetition frequency of 16 Hz. The calculations are for 
the temperature history at a position on the beam axis, 10|Lim beneath the surface of the tear layer; 
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Figure 2. Calculated temperature histories on the beam axis 10 pm 
beneath the tear surface at the threshold exposure for sequences of 32 
pulses, and (b) 128 pulses at 16 Hz. 
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thus, assuming a tear layer thickness of about 7 |xm, the temperature history is that just inside the 

anterior-most epithelial cells.Using the modified program we showed that the maximum 
temperature increase occurred -164 p,sec after the final pulse. The maximum temperature increases 
for all of the damage thresholds are listed in Table 1. 

The damage condition is described by an essentially constant maximum temperature rise 
and is consistent with a critical peak temperature damage model. For the exposures at 10 Hz, 
damage occurs at = 29.4±2.8 C (mean±SD), and for the exposures at 16 Hz, it occurs at 
= 32.0±2.0 C. TTie maximum temperature rises for the two exposure conditions are shown 
graphically in Figure 3. 
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Figure 3. The calculated maximum temperature rises on the beam axis 
10 pm beneath the tear surface for the damage threshold exposure. The 
lines show the mean values of for the two exposure conditions. 


In the Introduction we noted that light and electron micrographs of corneas exposed to a 
single 80 ns pulse just above the damage threshold show features consistent with both thermal and 
mechanical (e.g., acoustic) damage.^- 7 However, the fact that the calculated maximum 
temperature rises listed in Table 1 and plotted in Figure 3 are essentially constant for all of damage 
thresholds suggests that, at least for the multiple-pulse exposures, the damage mechanism may be 
predominately thermal. 

In order to test this suggestion, we performed a series of damage experiments on 
enucleated eyes cooled to room temperature (average 21 C) as described in Methods. The 
underlying hypothesis for this test is that if the critical temperature model is valid, then damage 
should occur for exposures that result in the same final critical temperature (not temperature 
increase). Thus in a cornea initially at room temperature, sufficient additional energy would have 
to be supplied first to raise the temperature to the in vivo temperature and then to the final d^age 
temperature. In a preliminary experiment using 8 pulses at 16 Hz we showed that corneas in the 
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cooled enucleated eyes did not incur damage at the same (slightly above threshold ) exposure that 
produced a lesion in the in vivo cornea. This finding suggested that the damage mechanism was at 
least partially thermal. We then proceeded to determine damage thresholds in cooled corneas for 
sequences of 8 and 32 pulses at 16 Hz. The threshold energy densities and the calculated 
maximum temperature rises for these experiments are listed in Table 2. Assuming that the 
temperature of the anterior surface of an in vivo cornea is 35 C, the “damage temperatures” from 
Table 1 for the 8 and 32 pulse 16 Hz thresholds are respectively 339.5 K and 341.1 K; whereas 
they are respectively 354 K and 345.8 K in the corneas cooled to 21 C before exposure. The 
additional energy required to produce a minimal lesion in the cooled corneas is therefore sufficient 
to raise the cornea temperature to a level slightly higher than that associated with damage in vivo. 
These results are highly suggestive that the multiple pulse damage is indeed predominately thermal. 
We speculate that the increased “damage temperatures” in the enucleated eye experiments may be 
due a slowing down of the processes leading to the observed damage endpoint as a result of the 
lower ambient temperature. (Recall that we determine damage 1/2 hour after exposure and the 
enucleated eye is maintained in BSS at 21 C during this time.) 

Table 2: Threshold Energy Densities and Calculated Maximum 
Temperature Rises for Enucleated Eyes at 21 C. 







8 

16 

393 

3.58 

59.8 

32 

16 

236 

3.58 

51.6 


*Calculated on the beam axis 10 pm beneath the anterior tear surface. 

We have submitted corneas with lesions to Dr. W. R. Green at the Wilmer Institute for 
histology. Corneas exposed slightly above the damage threshold for 8 pulses and 1024 pulses at 
16 Hz were submitted prior to the end of this reporting period; however no results were available. 
We plan to submit additional tissue exposed just above the damage threshold for 8 pulses, 128 
pulses, and 1024 pulses at 10 Hz and 128 pulses at 16 Hz. 

Conclusions 

Threshold damage to the corneal epithelium resulting from exposure to sequences of 80 ns 
pulses of CO 2 laser radiation is correlated by a power law of the form = CiV“ in which 
is the threshold energy density and N is the number of pulses in the sequence. The constant C and 
exponent a appear to differ slightly depending on the pulse repetition frequency. For sequences 
of pulses at 10 Hz, C = 291 mJ/cmVpulse and a = 0.162 and for sequences at 16 Hz, C = 300 

mJ/cm^/pulse and a = 0.194. Temperature calculations reveal that the maximum temperature 
increase on the beam axis 10 ^m beneath the anterior tear surface resulting from the different 
threshold exposures is essentially constant. For the exposures at 10 Hz, damage occurs at = 
29.4±2.8 C (mean+SD), and for the exposures at 16 Hz, it occurs at = 32.0±2.0 C. This 
result is consistent with a critical peak temperature damage model and suggests that, at least for the 
multiple-pulse exposures, the damage mechanism may be predominately thermal. D^age 
threshold measurements on corneas maintained at 21 C indicated that the damage mechanism is 
indeed predominately thermal. 
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APPENDIX 2 


Modified program for calculating temperature histories for sequences of very short pulses. 
Program allows for variable resolution following pulses. 


C PROGRAM CALCULATES TEMPERATURE-TIME HISTORIES FOR RADIATION 
INCIDENT 

C ON A SEMI-INFINITE ABSORBING SLAB 
C RADIATION CAN BE EITHER CW OR SEQUENCES OF PULSES 
C RADIATION AND CONVECTIVE HEAT TRANSFER ARE IGNORED 
C PROGRAM ADAPTED BY R. MCCALLY SEPTEMBER 25, 1997 FROM PROGRAM 
TCAL.F TO PROVIDE INCREASED TIME 
C RESOLUTION AFTER PULSES 
C 

C SO =IRRADIANCE(W/(CM*CM)) 

C KAPPA = THERMAL CONDUCTIVITY (CAL/(CM*DEGC)) 

C R = RADIAL DIST FROM BEAM AXIS (CM) 

C Z = DEPTH (CM) 

C RE = 1/E BEAM RADIUS (CM) 

C TAU = PULSE DURATION (SEC) 

C TAAFT = INCREASED RESOLUTION INTERVAL AFTER PULSES 
C C = HEAT CAPACITY (CAL/(GM*DEGC)) 

C GAM = INVERSE ABSORPTION LENGTH (1/CM) 

C DELT = 1/PRF (SEC) — DELT MUST BE >= TAU 
C NP = NUMBER OF PULSES 

C ND = NUMBER OF INTERVALS CALCULATED DURING EACH PULSE 
C NDAP = NUMBER OF POINTS IN THE INTERVAL TAAFT AFTER EACH PULSE 
C NPTS = NUMBER OF CALCULATED POINTS = NP*(ND+NDAP) -l- 1 (CALC IN 
PGM) 

NRUNS = NUMBER OF SETS OF INPUT DATA 

IRULE - CHOICE OF QUADRATURE RULE IN DQDAG. (INPUT) 

A GAUSS-KRONROD RULE IS USED WITH 
7 - 15 POINTS IF IRULE = 1 
10 - 21 POINTS IF IRULE = 2 
15 - 31 POINTS IF IRULE = 3 
20 - 41 POINTS IF IRULE = 4 
25-51 POINTS IF IRULE = 5 
30 - 61 POINTS IF IRULE = 6 
IRULE = 2 IS RECOMMENDED FOR MOST FUNCTIONS. 

IF THE FUNCTION HAS A PEAK SINGULARITY USE IRULE = 1. 

IF THE FUNCTION IS OSCILLATORY USE IRULE = 6. 

PROGRAM IRTEMP 
C 

IMPLICIT REAL*8(A-H,0-Z) 

REAL*8 KAPPA 

COMMON/BLOCKA/ ZO,ONE,TWO,TEN,HUN,SPLW,ZA,GTA,AZ,E2,R2, XO, 
TEGZ, G2FA 
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DIMENSION TIME(0:62000), Tl(0:62000), TEMP(0:62000),TEMPN(0:62000) 

C T1 ARE TIME POINTS DURING AND AFTER A SINGLE PULSE 
CHARACTER*! TAB 
EXTERNAL FF 
C 

INN = 7 ! DEFINES INPUT LOGICAL UNIT 

lOUT = 8 ' DEFINES OUTPUT LOGICAL UNIT FOR FINAL DATA FILE 

ISCREEN = 6 ! DEFINES LOGICAL UNIT FOR OUTPUT TO SCREEN 

ICOUNT = 0 'INITIALIZE COUNTER FOR TIMES INPUT FILE IS ACCESSED 


C 

C 


C 

c 


DEFINE CONSTANTS 
ZO = O.ODO 
ONE= l.ODO 
TWO = 2.0D0 
TEN = 1 O.ODO 
HUN = lOO.ODO 
PI = 4.0DO * DATAN( ONE) 

SPI = DSQRT( PI) 

TAB = CHAR(9) ! TAB IS ASC 9 


OPEN I/O FILES 

OPEN (UNIT=INN, FILE ='IRINHRES', STATUS='OLD',ACCESS='SEQUENTIAL, 
BLANK='NULL’) 

OPEN (UNIT=IOUT, FILE =TROUT.HIRES, 

STATUS='UNKNOWN',ACCESS='SEQUENTIAL’) 


C 

1000 CONTINUE 
C 

C GET INPUT DATA 


INPUT(SO,KAPPA,R,Z,RE,TAU,TAAFT,C,GAM,DELT,NP,ND,NDAP,RHO,NRUNS,I 

COUNT, 

1 INN,lOUT,ISCREEN, AERR,RERR,IRULE) 

C 

NT = ND+NDAP 


NPTS = NP*(NT+1) 
WRITE(ISCREEN,3)NPTS 
WRITE(IOUT,3)NPTS 

A2 = RHO*C/(4.0D0*KAPPA) 
TSIG = RE*RE 


R2 = R*R 

FRONT = SO*GAM/(8.36DO*RHO*C) 
ZA = Z**2 * A2 
AZ = DSQRT(ZA) 

TA = TWO * DSQRT( A2 ) 

GTA = GAM/T A 
XO = TWO*A2*Z/GAM 
TEGZ = TWO * DEXP(-GAM*Z) 
G2FA = GAM**2 / (4.D0*A2) 

W = A2 * TSIG 


C 

C CALCULATE TIME POINTS 
DO 30 I=0,ND 
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T1(I) = FTAU/ND 

30 CONTINUE 
C 

NDl =ND+ 1 
C 

DO 31 I=ND1,NT 

T1(I) = TAU + (I-ND)*TAAFT/NDAP 

31 CONTINUE 
C 

IF(NP.EQ.l)GO TO 38 
C 

DO 36J=1,NP 

JNT = (J-1)*(NT+1) 

DO 32 I=JNT,JNT+ND 
K = I - JNT 

TIME® = T1(K) + (J-1)*DELT 

32 CONTINUE 

DO 33 I=JNT+ND1,JNT+NT 
K = I - JNT 

TIME® = T1(K) + (J-1)*DELT 

33 CONTINUE 
36 CONTINUE 

C 

GO TO 40 
C 

38 DO 391=0,NT 

TIME® = T1(I) 

39 CONTINUE 
C 

40 TIME(NPTS) = NP*DELT 
C 

C BEGIN TEMPERATURE CALCULATION 
TEMP(O) = ZO 
C WRITE(ISCREEN,4) 

C WRITE(IOUT,5)TAB,TAB,TAB,TAB 
C 

DO 100I=1,NPTS 
TL = TIME(0) 

IF(TIME(I)-TAU .GT. TL) TL=TIME®-TAU 
TU = TIME(I) 

CALL DQDAG (FF, TL, TU, AERR, RERR, IRULE, PSI, ERREST) 

C DQDAG IS AN IMSL ROUTINE FOR EVALUATING THE INTEGRAL. IT IS A 
C GLOBALLY ADAPTIVE SCHEME THAT IS BASED ON GAUSS-KONRAD 

C RULES 

TEMP(I) = FRONT * PSI 

C WRITE(ISCREEN,6) I,PSI,TIME(I),TEMP(I),ERREST 
C WRITE(IOUT,7) I,TAB,PSI,TAB,TIME(I),TAB,TEMP(I),TAB,ERREST 
100 CONTINUE 
C 

DO 150 I=0,NPTS 

TEMPN(I) = TEMP® 

150 CONTINUE 
IF (NP.EQ.1) GO TO 375 
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C 

C HANDLE MULTIPLE PULSES 
C 

DO 300 J=2,NP 

JNT = (J-1)*(NT+1) 

DO 200 I=0,NPTS 

IF (I.LT.JNT) THEN 

TEMPN(I) = TEMPN(I) 

ELSE 

TEMPN(I) = TEMPN(I) + TEMP(I-JNT) 

END IF 

200 CONTINUE 

300 CONTINUE 
C 

C WRITE TO FINAL DATA FILE 

375 WRITE(IOUT,380)TAB,TAB ! WRITE HEADER FOR RESULTS 

DO 400 I=0,NPTS 

WRITE(IOUT,390) I,TAB,TIME(I),TAB,TEMPN(I) 

400 CONTINUE 

IF (ICOUNT .LT. NRUNS) GO TO 1000 
99 STOP 

3 FORMATC NPTS = ’,15,' TOTAL POINTS ') 

4 FORMAT(T3,T,T 11 ,'PSr,T25,’TIME',T36,'TEMPERATURE’,T52,'ERROR') 

5 FORMAT(T,Al,’PSr,Al,'TIME',Al,'TEMPERATURE’,Al,’ERROR') 

6 F0RMAT(I5,1PG14.6,1PG14.6,1PG14.6,1PG14.6) 

7 F0RMAT(I5,A1,1PG14.7,A1,1PG14.7,A1,1PG14.7,A1,1PG14.7) 

380 FORMAT('r,Al ,'TIME’,A1 ,'TEMPERATURE') 
390FORMAT(I5,A1,1PE14.8,A1,1PE14.8) 

END 

C 

C 

FUNCTION FF ( X ) 

C THIS DESCRIBES THE INTEGRAND 
IMPLICIT REAL*8(A-Z) 

COMMON/BLOCKA/ ZO,ONE,TWO,TEN,HUN,SPI,W,ZA,GTA,AZ,E2,R2,XO, 
ITEGZ, G2FA 
C 

IF(X .EQ. ZO) GO TO 50 
C 

X2 = DSQRT( X ) 

XW = ONE + XAV 
EE = -R2/XW 
El =ZO 

IF(DABS(EE) .LE. HUN) El = DEXP( EE ) 

C 

E2 = -ZA/X 
C 

Cl = GTA * X2 
C2 = AZ/X2 
ARGl = Cl + C2 
ARG2 = Cl - C2 
HI = HFUN( ARGl) 

H2 = HFUN( ARG2 ) 

C 
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IF(X .GT. XO) GO TO 10 

C... FOR X LESS THEN XO AND NOT ZERO 6/2/80 THE REVISION 
FO = El/XW 

G1 = TEGZ * DEXP(G2FA*X) 

G2=ZO 

IF(DABS(E2) .GT. HUN) GO TO 30 
HI = HFUN( DABS(ARGl)) 

H2 = HFUN( DABS(ARG2)) 

G2 = (H1-H2) 

C 

30 FF = FO*(Gl + G2) 

RETURN 

C 

10 CONTINUE 
FF = E1 * (HI +H2)/XW 
RETURN 

C... FOR X= 0 ONLY 
50 CONTINUE 
FF = TEGZ * DEXP(-R2) 

RETURN 

END 


FUNCTION HFUN ( ARG ) 

IMPLICIT REAL*8(A-Z) 

COMMON/BLOCKA/ ZO,ONE,TWO,TEN,HUN,SPLW,ZA,GTA,AZ,E2,R2, XO, 
ITEGZ, G2FA 

DATAHALF,ONEH/0.5D0, 1.50D0/ 

C 

Y = DABS(ARG) 

Y2 = Y*Y 
HFUN = ZO 
BB = Y2 + E2 

IF(ARG .GE. -13.0D0) GO TO 10 
Bl=ZO 

IF(DABS(BB) .LE. 174.0D0) B1 = TWO * DEXP( BB ) 

B2 = ZO 

IF(DABS(E2) .LT. 174.0D0) B2 = DEXP(E2) / 

* (SPF(Y+HALF/(Y+ONE/(Y+ONEH/(Y+TWO/Y))))) 

HFUN = B1 -B2 
RETURN 
10 CONTINUE 
IF(ARG .GE. ZO) GO TO 20 

IF(DABS(BB) .LT. 174.0D0) HFUN = DEXP(BB)*(TWO-DERFC(Y)) 

RETURN 
20 CONTINUE 
IF(ARG .GT. 13.0) GO TO 30 

IF(DABS(BB) .LT. 174.0D0) HFUN = DEXP(BB)*DERFC(ARG) 

RETURN 
30 CONTINUE 

IF(DABS(E2) .GT. 174.0D0) GO TO 40 
Z = ARG 

HFUN = SPP(Z+HALF/(Z+ONE/(Z+ONEH/(Z+TWO/Z)))) 

IF(HFUN .NE. ZO) HFUN = DEXP(E2)/HFUN 
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40 RETURN 
END 


FUNCTION DERFC(X) 

RETURNS THE COMPLIMENTARY ERROR FUNCTION 

CHECKS WITH MATHEMATIC A TO 15 SIGNIFICANT HGURES 
IN RANGE (.001 <= X <= 4.0) (MAXIMUM PRINTOUT FROM MATHEMATICA 
AND MAXIMUM RANGE CHECKED 3/18/96 RLM) 

IMPLICIT REAL*8 (A-H,0-Z) 

DATA Z0,HALF,ONE/0.0D0,0.5D0,1 .ODO/ 

IF(X.LT.Z0)THEN 

DERFC=ONE+GAMMP(HALF,X**2) 

ELSE 

DERFC=GAMMQ(HALF,X**2) 

ENDIF 

RETURN 

END 

FUNCTION GAMMP(A,X) 

IMPLICIT REAL*8 (A-H,0-Z) 

DATA Z0,HALF,ONE/0.0D0,0.5D0,1 .ODO/ 

IF(X.LT.ZO.OR.A.LE.ZO)PAUSE 
IF(X.LT.A+ONE)THEN 
CALL GSER(GAMSER,A,X,GLN) 

GAMMP=GAMSER 

ELSE 

CALL GCF(GAMMCF,A,X,GLN) 

GAMMP=ONE-GAMMCF 

ENDIF 

RETURN 

END 

FUNCTION GAMMQ(A,X) 

IMPLICIT REAL*8 (A-H,0-Z) 

DATA Z0,HALF,ONE/0.0D0,0.5D0,1 .ODO/ 

IF(X.LT.ZO.OR.A.LE.ZO)PAUSE 
IF(X.LT.A+ONE)THEN 
CALL GSER(GAMSER,A,X,GLN) 

GAMMQ=ONE-GAMSER 

ELSE 

CALL GCF(GAMMCF,A,X,GLN) 

GAMMQ=GAMMCF 

ENDIF 

RETURN 

END 


SUBROUTINE GSER(GAMSER,A,X,GLN) 
IMPLICIT REAL*8 (A-H,0-Z) 

DATA Z0,HALF,ONE/0.0D0,0.5D0,1 .ODO/ 
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PARAMETER (ITMAX=400,EPS=10D-17) 

GLN = 0.572364942924700IDO ! LN(SQRT(PI)) = GAMMA( 1/2) 

IF(X.LE.ZO)THEN 
IF(X.LT.ZO)PAUSE 
GAMSER=ZO 
RETURN 
END IF 
AP=A 

SUM=ONE/A 
DEL=SUM 
DO 11 N=1,ITMAX 
AP=AP+ONE 
DEL=DEL*X/AP 
SUM=SUM+DEL 

IF(DABS(DEL).LT.DABS(SUM)*EPS)GO TO 1 
11 CONTINUE 

PAUSE 'A TOO LARGE, ITMAX TOO SMALL 

I GAMSER=SUM*DEXP(-X+A*DLOG(X)-GLN) 

RETURN 

END 

C 

SUBROUTINE GCF(GAMMCF,A,X,GLN) 

IMPLICIT REAL*8 (A-H,0-Z) 

DATA Z0,HALF,ONE/0.0D0,0.5D0,l .ODO/ 

PARAMETER (ITMAX=400,EPS=:10D-17) 

GLN = 0.5723649429247001DO ! LN(SQRT(PI)) = GAMMA( 1/2) 

GOLD=ZO 
AO=ONE 
A1=X 
B0=Z0 
Bl=ONE 
FAC=ONE 
DO 11N=1,ITMAX 
AN=DFLOAT(N) 

ANA=AN-A 

A0=(A1+A0*ANA)*FAC 
B0=(B 1+B0*ANA)*FAC 
ANF=AN*FAC 
A1=X*A0+ANF*A1 
B1=X*B0+ANF*B1 
IF(ALNE.Z0)THEN 
FAC=ONE/Al 
G=B1*FAC 

IF(DABS((G-GOLD)/G).LT.EPS)GO TO 1 
GOLD=G 
ENDIF 

II CONTINUE 

PAUSE 'A TOO LARGE, ITMAX TOO SMALL’ 

1 GAMMCF=DEXP(-X+A*DLOG(X)-GLN)*G 
RETURN 
END 
C 

C- 

SUBROUTINE INPUT(SO,KAPPA,R,Z,RE,TAU,TAAFT,C,GAM,DELT,NP,ND, 
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1NDAP,RH0,NRUNS,IC0UNT,INN,I0UT,ISCREEN,AERR,RERR,IRULE) 

C- 

IMPLICIT REAL*8(A-H,0-Z) 

REAL*8 KAPPA 

C INTEGER*4 NPTS,NP,ND,NDAP,NRUNS,ICOUNT,INN,IOUT,ISCREEN,IRULE,I 
CHARACTER*9 DAY 
CHARACTER* 11 ADI 
CHARACTER*30 AIN2 

PURPOSE: READ INPUT FILE IRINPUT AND WRITE PARAMETERS. 

INPUT: INN,IOUT,ISCREEN 


OUTPUT: SO,KAPPA,R,Z,RE,TAU,TAAFT,C,GAM,DELT,NP,ND,NDAP, 
RHO,NRUNS,ICOUNT,AERR,RERR,IRULE. 

ICOUNT=ICOUNT+l ! INCREMENT FILE ACCESS COUNTER 

READ FILE INN. 

READ(INN,'(A30)',END=99,IOSTAT=I)AIN2 ! READ HEADER 

READ(INN,’(A30)’,END=99,IOSTAT=I)AIN2 ! READ — 

READ(INN,'(A11)',END=99,I0STAT=I)AD1 ! READ123 ETC 

READ(INN,'(A11,15)',END=99,I0STAT=I)AD1 ,NRUNS 
READ(INN,’(A11,1PG14.6)’,END=99,IOSTAT=I)AD1,SO 
READ(INN,'(A11,1PG14.6)’,END=99,I0STAT=I)AD1,TAU 
READ(INN,’(A11,1PG14.6)’,END=99J0STAT=I)AD1,TAAFT 
READ(INN,'(A 11 ,F 14.9)',END=99,IOSTAT=I) ADI ,RE 
READ(INN,'(A11 ,F 14.9)’,END=99,IOSTAT=I)AD 1 ,Z 
READ(INN,'(A 11 ,F14.9)’,END=99,IOSTAT=I)AD 1 ,R 
READ(INN,’( A11,15)’,END=99,IOSTAT=I) AD 1 ,NP 
READ(INN,'(A11,1PG14.6)',END=99,I0STAT=I)AD1 ,DELT 
READ(INN,'( A11,15)',END=99,IOSTAT=I) AD 1 ,ND 
READ(INN,'(A11,15)',END=99,IOSTAT=I)AD 1 ,NDAP 
READ(INN,’(All,F14.9)',END=99,IOSTAT=I)ADl,GAM 
RE AD(INN,’( A11 ,F 14.9)',END=99,IOST AT=I)AD 1 ,C 
READ(INN,’(All,F14.9)',END=99,IOSTAT=I) ADI,KAPPA 
READ(INN,'(All,F14.9)',END=99,IOSTAT=I)ADl,RHO 

.INTEGRATION PARAMETERS: 

READ(INN,'( A11 ,F 14.9)’,END=99,IOSTAT=I)AD 1 ,AERR 
READ(INN,'(A 11, IPDl 4.9)',END=99,IOSTAT=I)ADl ,RERR 
READ(INN,'(A11,15)’,END=99,I0STAT=I)AD1 ,IRULE 

WRITE INPUT TO FILE IROUT ( DISK DATA FILE) 

WRITE(IOUT,100) 

WRITE(IOUT,105) 

CALLDATE(DAY) 

WRITE(IOUT,l 10)DAY 
AIN2 = ■ INPUT PARAMETERS' 

WRITECIOUT, 110)AIN2 
WRITECIOUT, 150)NRUNS 
WRITE(IOUT,151)SO 
WRITECIOUT, 152)TAU 
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WRITE(IOlJT, 167)TAAFT 
WRITE(IOUT,153)RE 
WRITE(IOUT,154)Z 
WRITE(IOUT,166)R 
WRITE(IOUT,155)NP 
WRrrE(IOUT, 156)DELT 
WRITE(IOUT,157)ND 
WRITE(IOUT,158)NDAP 
WRITE(IOUT, 159)GAM 
WRITE(IOUT,160)C 
WRITE(IOUT, 161 )KAPP A 
WRITE(IOUT, 162)RHO 
AIN2 = ' INTEGRATION PARAMETERS: ' 
WRITE(IOUT, 110)AIN2 
WRrrE(IOUT, 163)IRULE 
WRrrE(IOUT, 164) AERR 
WRITE(IOUT, 165)RERR 
C 

C WRITE INPUT TO SCREEN 
C 

WRITE(ISCREEN, 100) 

WRITE(ISCREEN,110)DAY 
AIN2 = ’ INPUT PARAMETERS' 
WRITE(ISCREEN,110)AIN2 
WRITEdSCREEN, 150)NRUNS 
WRITE(ISCREEN, 151 )S0 
WRITE(ISCREEN,152)TAU 
WRITE(ISCREEN,167)TAAFr 
WRITEdSCREEN,153)RE 
WRITEdSCREEN,154)Z 
WRITE(ISCREEN, 166)R 
WRITE(ISCREEN, 155)NP 
WRITE(ISCREEN, 156)DELT 
WRITEdSCREEN, 157)ND 
WRrrE(ISCREEN,158)NDAP 
WRITEdSCREEN, 159)GAM 
WRITEdSCREEN, 160)C 
WRITEdSCREEN, 161 )KAPP A 
WRITEdSCREEN, 162)RHO 
AIN2 =' INTEGRATION PARAMETER: ' 
WRITE(ISCREEN,110)AIN2 
WRITEdSCREEN, 163)IRULE 
WRITE(ISCREEN, 164)AERR 
WRITEdSCREEN, 165)RERR 

99 CONTINUE 
C 

C FINAL WRITE STATEMENTS 
C 

WRITEdOUT,120) 

WRITEdSCREEN, 120) 

RETURN 

C 

100 FORMAT(/' —WE ARE IN IRINHRES — '/) 
105 FORMAT^ — FILE IROUT.HIRES— '/) 


23 



DAMD17-96-C-6005 


no FORMAT(A35) 

120 FORMAT(/' —LEAVING IRINHRES— ’/) 

130 FORMATC ... ') 

150 FORMATC NRUNS =',15) 

151 FORMATC SO =’,1PG14.5,'WATTS/CM**2') 

152 FORMATC TAU = ',1PG14.5,' SEC) 

153 FORMATC RE = ',F14.5,’ CM’) 

154 FORMATC Z = ’,F14.5,' CM') 

155 FORMATC NP =',15,' PULSES') 

156 FORMATC DELT = ',1PG14.5,' SEC) 

157 FORMATC ND =',15,' POINTS/PULSE’) 

158 FORMATC NDAP =',15,' POINTS AFTER PULSE') 

159 FORMATC GAM = ',F14.5,' 1/CM') 

160 FORMATC C = ’,F14.5,' CAL/(GM*DEGC)') 

161 FORMATC KAPPA = ',F14.5,' CAL/(CM**2*DEGC)') 

162 FORMATC RHO = ',F14.5,' GM/CM**2') 

163 FORMATC IRULE =',15) 

164 FORMATC AERR = ',F14.5) 

165 FORMATC RERR = ',1PD14.5) 

166 FORMATC R = ',F14.5,' CM ) 

167 FORMATC TAAFT =’,1PG14.5,'SEC) 

END 
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